#include <stdio.h>
#include <math.h>

/* 返回值在区间： [-1,1]          */
/* 如返回-10，则证明输入参数无效    */
#define delta 0.0001f
double calculate_corss_correlation(double *s1, double *s2,int n)
{
    double sum_s12 = 0.0;
    double sum_s1  = 0.0;
    double sum_s2  = 0.0;
    double sum_s1s1 = 0.0; //s1^2
    double sum_s2s2 = 0.0; //s2^2
    double pxy = 0.0;
    double temp1 = 0.0;
    double temp2 = 0.0;
    
    if( s1==NULL || s2==NULL || n<=0)
      return -10;
    
    for(int i=0;i<n;i++)
    {
         sum_s12  += s1[i]*s2[i];
         sum_s1   += s1[i];
         sum_s2   += s2[i];
         sum_s1s1 += s1[i]*s1[i]; 
         sum_s2s2 += s2[i]*s2[i]; 
    }
    
    temp1 = n*sum_s1s1-sum_s1*sum_s1;
    temp2 = n*sum_s2s2-sum_s2*sum_s2;
    
    /* 分母不可为0 */
    if( (temp1>-delta && temp1<delta) || 
        (temp2>-delta && temp2<delta) ||
        (temp1*temp2<=0) )
    {
        return -10;  
    }       
    
    pxy = (n*sum_s12-sum_s1*sum_s2)/sqrt(temp1*temp2);
    
    return pxy;
}
double s1[30] = {
0.309016989,0.587785244,0.809016985,0.95105651,1,0.951056526,
0.809017016,0.587785287,0.30901704,5.35898E-08,0,0,
0,0,0,0,0,0,
0,0,0,0,0,0,
0,0,0,0,0,0
};
double s2[30] = {
0.343282816,0.686491368,0.874624132,0.99459642,1.008448609,
1.014252458,0.884609221,0.677632906,0.378334666,0.077878732,
0.050711886,0.066417083,0.088759401,0.005440732,0.04225661,
0.035349939,0.0631196,0.007566056,0.053183895,0.073143706,
0.080285063,0.030110227,0.044781145,0.01875573,0.08373928,
0.04550342,0.038880858,0.040611891,0.046116826,0.087670453
};

int main(void)
{
    double pxy;
    double s3[30];

    pxy = calculate_corss_correlation(s1,s2,30);
    printf("pxy of s1 and s2:%f\n",pxy);  
    
    pxy = calculate_corss_correlation(s1,s1,30);
    printf("pxy of s1 and s1:%f\n",pxy);  

    for(int i=0;i<30;i++)
    {
        s3[i] = -1*s1[i];
    }
    pxy = calculate_corss_correlation(s1,s3,30);
    printf("pxy of s1 and s3:%f\n",pxy);  
    
    return 0;
}
